################################################
# Analysis for:
#
# Haglin, Kathryn, Soren Jordan, Alison Higgins Merrill, and Joseph Daniel Ura. 
#  ``Reevaluating Ideological Asymmetries in Specific Support for the Supreme Court.''
#  Research & Politics
#
# Required files: HJMU-RP-2025.dta
#
################################################

library(tidyverse)
library(haven)
library(dynamac)
library(urca)
library(systemfit)
library(viridis)
library(tseries)
library(broom)
library(tidyverse)
library(car)
library(lmtest)
library(colorspace)
library(plotly)
library(lattice)
library(ggpubr)


data <- read_dta("/.../HJMU-RP-2025.dta") 

data
tail(data) # be sure it includes 2024




##################################################################################################################
# Figure 1. Overall trends + Approval gap
##################################################################################################################

data.viz.all <- data %>% 
	select(year, tooliberalmacro, tooconservativemacro, aboutrightmacro) %>%
	pivot_longer(cols = tooliberalmacro:aboutrightmacro) 

data.viz.all$Evaluation <- NA
data.viz.all$Evaluation[data.viz.all$name == "tooliberalmacro"] <- "Too Liberal"
data.viz.all$Evaluation[data.viz.all$name == "tooconservativemacro"] <- "Too Conservative"
data.viz.all$Evaluation[data.viz.all$name == "aboutrightmacro"] <- "About Right"
data.viz.all$Evaluation <- factor(data.viz.all$Evaluation, levels = c("Too Conservative", "Too Liberal", "About Right"))


evals <- ggplot(data.viz.all) +
  geom_line(aes(x = year, y = value, color = Evaluation, linetype = Evaluation), linewidth = 1.5, alpha = 0.8) +
  # geom_line(data = as.data.frame(data), aes(x = year, y = approvescmacro), linewidth = 2.5) +
  scale_color_manual(values = c("darkred", "darkblue", "grey50")) +
  # geom_text(x = 2017, y = 60, label = "Overall Approval", size = 6) + 
  scale_x_continuous(breaks = seq(2000, 2024, 4)) +
  scale_y_continuous(limits = c(10, 55)) +
  annotate("text", x = 2002.5, y = 53, label = "About Right", size = 6) +
  annotate("text", x = 2022.5, y = 11, label = "Too Liberal", size = 6) +
  annotate("text", x = 2002.5, y = 16, label = "Too Conservative", size = 6) +
  labs(x = "Year", y = "Gallup Topline Percentage") +
  theme_bw() +
  theme(legend.position = "none") +
  theme(text = element_text(size = 20), 
  	# axis.text.x = element_text(angle = 45, hjust = 1),
  	strip.text.x = element_text(size = 14), strip.text.y = element_text(size = 14))
# dev.off()




data$upper <- data$lower <- NA
data$lower[data$year > 2019] <- data$res2018[data$year > 2019]
data$upper[data$year > 2019] <- 0

gap <- ggplot(data, aes(x = year, y = res2018)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey", alpha = 0.5) +
  scale_y_continuous(limits = c(-20, 10)) +
  geom_line(linewidth = 1.5) +
  scale_x_continuous(breaks = seq(2000, 2024, 4)) +
  geom_hline(yintercept = 0, lwd = 1.5) +
  labs(x = "Year", y = "Approval Gap (Approval - Predicted)") +
  annotate("text", x = 2021.85, y = 6.5, label = "Out of Sample\nResiduals\n(Post-2018)", size = 6) +
  annotate("segment", x = 2018.5, xend = 2018.5, y = -10, yend = 10, linewidth = 1.25, alpha = 0.5) + 
  theme_bw() +
  theme(text = element_text(size = 20), 
  	# axis.text.x = element_text(angle = 45, hjust = 1),
  	strip.text.x = element_text(size = 14), strip.text.y = element_text(size = 14))

# pdf("/.../fig1-evalsgap.pdf", width = 10, height = 10)
ggarrange(gap, evals, nrow = 2)
dev.off()








##################################################################################################################
# Figure 2. SUR over time
##################################################################################################################
# Straight replicate PRQ
data$l.approvescmacro <- lshift(data$approvescmacro, 1)

sur.data <- data.frame(Year = 2015:2024,
						tl.est = NA,
						tl.low = NA,
						tl.high = NA,
						tc.est = NA,
						tc.low = NA,
						tc.high = NA)	
						
eqn1 <- approvescmacro ~ approvecongmacro + dempres + tooliberalmacro + l.approvescmacro
eqn2 <- approvescmacro ~ approvecongmacro + dempres + tooconservativemacro + l.approvescmacro
eqn3 <- approvescmacro ~ approvecongmacro + dempres + tooliberalmacro + tooconservativemacro + l.approvescmacro


for(i in 1:nrow(sur.data)) {
	fitsur <- systemfit(list(eqn1, eqn2, eqn3), data = data[data$year <= sur.data$Year[i],], method = "SUR")
	sur.data$tl.est[i] <- tidy(fitsur)$estimate[tidy(fitsur)$term == "eq3_tooliberalmacro"]
	sur.data$tl.low[i] <- tidy(fitsur)$conf.low[tidy(fitsur)$term == "eq3_tooliberalmacro"]
	sur.data$tl.high[i] <- tidy(fitsur)$conf.high[tidy(fitsur)$term == "eq3_tooliberalmacro"]
	sur.data$tc.est[i] <- tidy(fitsur)$estimate[tidy(fitsur)$term == "eq3_tooconservativemacro"]
	sur.data$tc.low[i] <- tidy(fitsur)$conf.low[tidy(fitsur)$term == "eq3_tooconservativemacro"]
	sur.data$tc.high[i] <- tidy(fitsur)$conf.high[tidy(fitsur)$term == "eq3_tooconservativemacro"]
}

						
# Poor version of pivot_longer
sur.data.long <- data.frame(Year = sur.data$Year,
					Estimate = c(sur.data$tl.est, sur.data$tc.est),
					Lower = c(sur.data$tl.low, sur.data$tc.low),
					Upper = c(sur.data$tl.high, sur.data$tc.high),
					Evaluation = c(rep("Too Liberal", length(sur.data$Year)), 
						rep("Too Conservative", length(sur.data$Year)))
					)

sur.data.long$Evaluation_reorder <- factor(sur.data.long$Evaluation, levels = c("Too Liberal", "Too Conservative"))

# pdf("/.../fig2-sur-separate.pdf", width = 10, height = 7.5)
ggplot(sur.data.long, aes(x = Year, y = Estimate, color = Evaluation, shape = Evaluation)) +
	facet_grid(vars(Evaluation_reorder)) +
	geom_point(size = 4, position = position_dodge(width=0.4)) + 
	geom_errorbar(aes(ymin = Lower, ymax = Upper), 
		width = 0, linewidth = 1.5, position = position_dodge(width=0.3)) +
	scale_color_manual(values = c("darkred", "darkblue")) +
	scale_x_continuous(breaks = seq(2000, 2024, 4)) +
	theme_bw() + 
	labs(y = "Coefficient Estimate", x = "SUR Model from 2000 - Year") +
	theme(legend.position = "none") +
	theme(text = element_text(size = 20), 
		strip.text.x = element_text(size = 14), 
		strip.text.y = element_text(size = 14))
dev.off()









##################################################################################################################
# Supplemental Materials
##################################################################################################################


## For the SM: Tabular results for the SM. table of number of coefficients*2 rows x years columns
# Fit one SUR to obtain number of coefficients
fitsur <- systemfit(list(eqn1, eqn2, eqn3), data = data)
the.terms <- length(rep(tidy(fitsur)$term, each = 2))
the.years <- 2015:2024

the.surs <- matrix(rep(NA, the.terms*(length(the.years)+1)), nrow = the.terms)

the.surs[,1] <- rep(tidy(fitsur)$term, each = 2)

for(i in 1:length(the.years)) {
	fitsur <- systemfit(list(eqn1, eqn2, eqn3), data = data[data$year <= the.years[i],], method = "SUR")	
	the.surs[((1:(nrow(the.surs)/2))*2)-1, (i+1)] <- round(tidy(fitsur)$estimate, digits = 3)
	the.surs[((1:(nrow(the.surs)/2))*2), (i+1)] <- paste0("(", round(tidy(fitsur)$std.error, digits = 3), ")")
}

colnames(the.surs) <- c("Coef", the.years)
the.surs # prettify in latex; Table SM 2

# manually add the R2 since it would be annoying to break apart the equation results in the loop
for(i in 1:length(the.years)) {
	fitsur <- systemfit(list(eqn1, eqn2, eqn3), data = data[data$year <= the.years[i],], method = "SUR")	
	print(c(the.years[i], summary(fitsur)$eq[[1]]$adj.r.squared, summary(fitsur)$eq[[2]]$adj.r.squared, summary(fitsur)$eq[[3]]$adj.r.squared))
}












# Is the change from 2020 to 2021 in TC statistically significant?

tc.2021.est <- -1.180
tc.2020.est <- -0.384

tc.2021.se <- 0.204
tc.2020.se <- 0.256

diff.2021 <- tc.2021.est - tc.2020.est
diff.se.2021 <- sqrt(tc.2020.se^2 + tc.2021.se^2)
diff.2021/diff.se.2021






